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ACCURATE  EIGENVALUES  OF  A  SYMMETRIC 


TRI-DIAGONAL  MATRIX 


By 

*/ 

W.  Ka.harf” 


ABSTRACT 

Having  established  tight  bounds  for  the  quotient  of  two  different 
lub-norms  of  the  same  tri-diagonal  matrix  J  ,  the  author  observes  that 
these  bounds  could  be  of  use  in  an  error -analysis  provided  a  suitable 
algorithm  were  found,  Such  an  algorithm  is  exhibted,  and  its  errors  are 
thoroughly  accounted  for,  including  the  effects  of  scaling,  over/ under¬ 
flow  and  roundoff.  A  typical  result  is  that,  on  a  computer  using  rounded 
floating  point  binary  arithmetic,  the  biggest  eigenvalue  of  J  can  be 
computed  easily  to  within  2,5  units  in  its  last  place,  and  the  smaller 
eigenvalues  will  suffer  absolute  errors  which  are  no  larger.  These 
results  are  somewhat  stronger  than  had  been  known  before, 
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Two  Questions : 


The  following  questions  are  connected  with  certain  error-analyses  of 
the  computed  eigenvalues  of  the  symmetric  tri-diagonal  NXN  matrix 


Our  notation  is  very  much  like  Householder's  (1964);  we  write 

lutg(A)  s  (max.  eigenvalue  of  AHA)1//2  and  lubg(A)  »  max^  | A^ ^  |  =  lulOp ( | A | )  , 

where  |a|^  s  |a^  |  .  The  questions  are 

1:  What  bounds  can  be  found  for  lubg(j)/lubg(j)  ? 

2:  What  bounds  can  be  found  for  lubg(j)/lubg( | j| )  ? 

We  shall  see  that  the  answers  are  respectively 

Is  ~\A<  lubs(J)/lubE(J)  <  lubs(j)/lubg(|j|)  . 


o  • 

C:  i 


~\ff: <  iubs(J)/lubs(  | j|  )  <  i  . 


The  only  new  results  here  are  the  lower  bounds  and  ~\j^ >  the  other 

inequalities  are  well  known  and  will  not  be  proved  here.  (For  proofs  see 
Householder's  book,  §§2.2  to  2.4,  with  which  the  reader  must  be  assumed 
to  have  extensive  acquaintance.)  Part  of  the  interest  in  the  constants 
~  and  ~\/“  arises  because  they  are  best  possible,  and  much  larger 


1 


than  the  lower  bound  ”\/1/N  which  would  be  required  if  J  were  replaced 
by  an  arbitrary  (symmetric)  matrix  in  the  inequalities  above.  A  brief 
survey  of  such  more  general  (and  therefore  weaker)  bounds  is  given  by 
Mrs.  B.  J.  Stone  (1962). 


Proof  of  1; 

The  results  which  we  wish  to  prove  are  insensitive  to  diagonal  simi¬ 
larity  transformations  and  to  the  replacement  of  J  by  -J  .  Therefore 
we  may  assume  without  loss  of  generality  that  all  b^  >  0  .  We  shall  write 

bo  3  bN  a  0  and 
ri  “  bi  +  bi-l  ' 

Hence 


lubE(j)  =  lubE(|j|)  »  maxi(|a1|  +  r^ 


for  some  k  defined  (perhaps  not  uniquely)  by  the  last  equation.  No 
generality  is  lost  by  assuming  that  a^  >  0  ,  so 


ak  +  rk  “  >  laj_l  +  ri  f°r  a11  i  • 

Now,  lubg(j)  is  the  largest  of  the  magnitudes  of  the  eigenvalues 
of  J  ,  so  the  minimax  characterization  of  those  eigenvalues  (see 
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Householder's  book,  §  3»3*l)  implies  that 


lubs(j)  >  lubs(K) 


for  any  principal  submatrix  K  of  J  .  It  is  particularly  convenient 
here  to  take 


K  2 


a 


k-1 

Vl 

0 


b, 


k-1 


b. 


0 


k+1 


A  related  matrix  K  is  obtained  by  reflecting  K  in  its  skew-diagonal; 


K  a 


\ 


a 


k+1 


b, 


0 


0 


k-1 


k-1  k-1 


This  reflection  changes  no  eigenvalue,  so 


lub0(K)  =  lubq(K)  >  —lub0 (K  +  £) 

o  o  —  d  o 


Consequently 


lubs(J)/lubE(j)  >  lubg(X) 

whe  re 

X  S  |(K  +  k)/(S]{  +  rk) 
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It  is  convenient  now  to  define 


x  ■  iV  (ak +  rk)  8nd 
y  '  <Vi  +  Vi^  (2*k  +  rk)  • 


Obviously  0  <  x  <  ^  .  Also  -1  <  y  <  1  because 


Vi  +  Vi>  *  IvJ  + 


^  <*k  +  rk  -  rk-x)  +  ("k  +  rk  •  Vi> 

<  (\  +  \  -  bk-1)  +  (ak  +  rk  -  bk) 


*  2alt  +  rk  * 


The  matrix  X  can  be  expressed  simply  in  terms  of  x  and  y  thus ; 


y-xy  x  0 

x  l-2x  x 


y-xy 


A  further  simplification  is  achieved  by  the  use  of  the  orthogonal  matrix 


0 

1 


V” 


vr 

0 

VF 


4 


and  the  2X2  matrix 


l-2x  V2x 


-X  y-xy 


which  are  connected  to  X  by  the  eigenvalue -preserving  similarity  trans¬ 
formation 


T 

QXQ  » 


z  o 

O  y-xy 


Since  y-xy  separates  the  two  eigenvalues  of  Z  ,  these  two  eigenvalues 
must  be  the  algebraically  greatest  and  least  eigenvalues  of  X  .  Therefore 
our  progress  so  far  can  be  summarized  by  the  inequality 


lub-(J)/luVi  (J)  >  lub.(Z)  , 


and  our  result  no.  1  will  be  proved  when  we  have  shown  that 


lubg(Z)  >  -\/i  . 


This  last  Inequality  is  obtained  below  from  a  demonstration  that 


7  =  min.lub  (Z)  over  (0  <  x  <  \  and  -1  <  y  <  1)  . 


Let  the  eigenvalues  of  Z  be  regarded  now  as  functions  of  y  for  a 
fixed  x  .  They  are  both  monotonic  non-decreasing  functions  of  y  because 


some  positive  multiple 


any  increase  in  y  is  tantamount  to  adding  to  Z, 
of  the  positive  semi~def in ite  matrix 

1  0  0 

I  0  1-x 


The  value 

yo  £  -  2x )/  ( 1  -  x) 


satisfies  -1  <  y  <  0  j  and  when  y  *  yQ  the  eigenvalues  of  Z  are 

just  +z  and  -z  ,  where 
o  o 

z  o  £  "W6x2  -  4x  4  1  , 

The  values  yQ  and  z^  are  significant  because  for  any  other  value  of  y 
the  matrix  Z  has  either  a  positive  eigenvalue  >  z 

or  a  negative  eigenvalue  <  •  z 

In  other  words, 


z  =  min  lub0 ( Z )  over  -1  <  y  <  .1. 


for  any  fixed  x  in  0  <  x  <  .  And  z * s  minimum  value 

when  x  -  -r  .. 

^  ,  rp 

The  foregoing  proof*  that  lub0  ( J )/  lu'b  ( J )  >  ~\l  -  also 

o  n  "  y  ; 

example 


iu  taken 


points  to  an 
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1 


0 


J  = 


1 

1 


with  lubg(J)  «  and  lubE(J)  «  3  j 
cannot  be  increased. 


therefore  the  lower  bound 


Proof  of  2: 

We  wish  to  show  that  lubg( J)/lubg(  |  j| )  > 
assume  without  loss  of  generality  that  all  b^ 
begin  with  some  definitions.  First  let 


As  before,  we 
It  is  convenient  to 


n  *  lubg(|j|) 

«  max.  eigenvalue  of  |j| 


Second,  define 


M  B  i(|j|  .  J)  . 

Ct* 


Evidently  M  is  a  non-negative  diagonal  matrix  whose  positive  elements 
are  just  the  positive  diagonal  elements  of  -J  .  For  the  sake  of  symmetry 
we  should  like  to  have  a  similar  definition  for  the  non-negative  diagonal 
matrix  P  whose  positive  elements  are  just  the  positive  diagonal  elements 
of  +J  .  Such  a  definition  is  provided  in  stages  as  follows,  We  define 

E  o  diag(-l,  +1,  -1, ...,  (-1)N)  , 
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and  use  it  in  an  eigenvalue-preserving  similarity  transformation  to  define 
the  matrix 

J  »  -EJE 

whose  eigenvalues  and  diagonal  elements  are  just  the  negatives  of  those 
of  J  .  But  J  has  the  same  off-diagonal  elements  as  J  and  |j|  . 
Therefore  the  matrix 

P  ■  |(|j|  -  J) 

is  defined  in  much  the  same  way  as  was  M  .  Note  that  PM  =  0  .  Finally, 
2  *"'2 

because  J  and  J  have  the  same  eigenvalues, 

\  ■  lub  (J)  =  lubs(J)  . 

Now  we  may  proceed  to  demonstrate  that  k/\i  > 

According  to  the  theory  of  non-negative  matrices  outlined  in  §2.4 
of  Householder '  s  book,  there  must  exist  some  non-negative  vector  v.  such 
that 

i  i  T 

|  J | y  (iv  >  0  and  v  v  =  1 

2  T  2  T 

Since  (lub^(J))  =  nx.x  J  x  over  x  x  -  1  , 

>  /A  =  VT(  I  jj  -  8M)2v 

2  l,  T„  ,  i,  T.,2 
=  (i  “  4|av  My  +  4v  M  v 
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Similarly, 


>  v^v  =  VT(  I  j|  -  2P)2V 

2  T  T  ? 

=  H  -  4|iv  Pv  +  4v  P  v 


Adding  and  using  the  fact  that  PM  =  0  yields 

2\2  >  2^2  -  WT(M  +  P)(til  -  M  -  P)v  . 
But  M  +  P  =  diag(|a^|)  ,  and 

|a±l  (m-  -  laj  )  <  • 

Therefore 

4v^(M  +  P) (pj  -  M  -  P)v  <  ji2v^v  =  \i~ 

and  so 

~2  .  n  2  2  2 
2a.  >  2\x  -  p,  =  p 


as  desired.  Result  no.  2  is  proved. 

This  proof  points  less  directly  than  did  the  proof  of  result  no.  1 
to  an  example  J  for  which  the  second  bound  is  achieved,  i.e.  for  which 

lubg(j)  « 
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In  fact,  the  foregoing  proof  was  motivated  by  a  foreknowledge  of  the  fol¬ 
lowing  example. 

Let  ai  =  |(-l)ix  for  1  <  i  <  N  ,  and  b.  =  |  for  1  <  i  <  N  . 

The  value  of  x  will  be  chosen  later  to  be  the  same  as  p  defined  above, 
but  first  we  observe  that  now 


P  +  M  =  |xl  ,  P  -  M  -  -jxE 


2  ' 


and 


c  *  |j|  -  4a  =  j  -  hts 


is  an  NXN  matrix  with  zero  on  the  diagonal  and  —  on  the  subdiagonal  and 
superdiagonal.  The  eigenvalues  of  C  are  well-known;  they  are  just  the 
numbers 


7n  s  cos  rut /  (H  +  l)  for  n  =  1,2,...,  N  . 


(See  Householder's  book,  p.  34  ex.  50.  His  matrix  J  is  defined  on  p.  2 
and  differs  from  ours.  His  K  =  2C  . )  In  particular,  since 


J|  *  C  +  gxl  , 


p  =  7i  + 


12* 


2 


Next  let  us  compute  the  largest  eigenvalue  \  of  J  The  computai ion 
is  considerably  shortened  by  Jim  Varah's  observation  that 


9  _2  1  2„ 

J  =  C  +  jjx  I  . 
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Therefore 


2.12 


o 

*  =  h  +  p 


The  ratio  \/n  takes  on  its  minimum  value 

T  fT 

2  when  x  «  \i  ~  ?.y ^  .  This  example  shows  that  the  lower  bound  y- 
cannot  be  increased. 


Application: 

Let  the  eigenvalues  k^  of  J  be  ordered  thus: 


\  <  x2  -  *  *  ’  -  Nl-l  -  S  5 


and  suppose  the  eigenvalues  (\^  +  bk i)  of  (J  +  6J)  are  ordered  simi¬ 
larly;  +  bk^  <  •  Here  the  matrix  5J  is  a  perturbation 

attributed,  possibly,  to  rounding  errors  in  a  numerical  calculation.  We 
shall  assume  that  5J  is  tri-diagonal  with  elements  bounded  by,  say, 


|&ail<a|ai|  and  |e>b± |  <  3 |b5 1 


where  a  and  p  are  small  positive  constants.  Given  a  and  p  ,  how 
big  can  be  ? 

The  easiest  bound  for  uses  the  fact  that,  if  the  eigenvalues  of 

5  J  a  re 

bl  1  &2  -  ' '  *  -  5N  ' 

then  6^  <  <  6^  . 
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A  proof  of  this  relation  can  be  found  in  Householder’s  book  (1964)  p.  79° 
Consequently 


|5\i|  <  lubg(6J)  <  lubg ( | S J | )  <  lubE(6J)  . 

In  particular,  if  a  =  £  then  | 5 J |  <  a| j|  elementwise,  so 

|8\J  <  yia  lut^(j)  =  y23  maxj|\^| 

for  all  l  by  virtue  of  the  inequality  no.  2.  More  generally,  inequality 
no.  2  can  be  extended  without  any  difficulty  to  the  case  that  (X  /  p  and 
yields  the  bound 


|a\d|  <  lubg( | 6j| )  <\{?  +  p2  lubg(  J )  . 

Though  pessimistic,  these  bounds  are  slightly  stronger  than  the  best  bounds 
available  in  terms  cf  lubg(J)  .  But  are  there  any  practical  circumstances 
where  such  bounds  may  be  of  use?  They  rely  upon  the  inequalities 

|Sa^|  <  a|a^|  and  1 6b  J  <  p|b.J  , 

whereas  the  typical  rounding  error  analyses  of  the  past  have  contained 
weaker  constraints  like 

1 6a.  I  <  a  lube(j)  and  |&b.|  <  plubJj) 

1  “  b  X  —  b 
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(cf.  Wilkinson's  book  (1965)  p.  304).  Thus  we  are  faced  with  the  following 
problem: 

Given  a  set  of  error-bounds,  find  a  numerical  algorithm  to  which 
they  are  applicable. 

This  problem  has  an  elegant  solution  which  is  described  below. 


The  Algorithm: 


We  shall  now  exhibit  and  completely  error- analyze  a  simple  and  effec¬ 
tive  method  for  computing  any  eigenvalue  of  J  .  The  basic  method 

was  first  put  forth  in  Dr.  Boris  Davison’s  numerical  analysis  lectures  at 
the  University  of  Toronto  in  1959>  and  begins  with 

Sylvester’s  Law  of  Inertia.: 

T 

Suppose  A  =  A  is  symmetric 

1  is  non-singular,  and 
DsL  bA(L  is  diagonal. 

Then  the  numbers  of  positive,  negative  and  zero 
diagonal  elements  of  D  are  the  same  respectively 
as  thv.  numbers  of  positive,  negative  and  zero  eigen¬ 
values  of  A  .  4 

A  proof  may  be  found  in  any  standard  text  on  matrices;  e.g.  in  Gantmacher 
(1959)  vol.  I  p.  297.  We  shall  apply  this  Law  to  the  triangular  factori¬ 
zation  of 


J  -  xl  -  LU  -  LDLT 


into  triangular  bi-diogonal  matrices  L  and  U  obtained  by  Gaussian 
elimination  without  pivotal  interchanges.,  It  is  unnecessary  to  compute 
any  but  the  diagonal  elements  u  of  U  .  They  are  obtained  from  the 
simple  recurrence 

u^  =  -  x  and 

p 

Un  =  an  "  X  "  bn-l//Vl  for  n  =  2>5'***>  N  • 
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but  such 


This  recurrence  breaks  down  if  and  only  if  some  value  u^  =  0  , 
a  thing  can  happen  only  if  x  takes  on  one  of  at  most  ~N(N  +  l)  excep¬ 
tional  values.  Indeed,  it  is  easy  to  see  that 

un  =  un(x)  *  <pn(x)/Vl(x) 

where  cpn(x)  the  characteristic  polynomial  of  the  first  n*n  principal 
submatrix  of  J  .  In  particular, 

<PN00  =  det(J  -  xl)  . 

Consequently,  the  recurrence  can  break  down  only  if  x  coincides  with 
one  of  the  eigenvalues  of  one  of  the  leading  principal  Liubmatrices  of  J  . 
Let  us  postpone  the  discussion  of  these  exceptional  values  of  x  ;  suppose 
for  now  that  the  recurrence  is  successful,  and  compute 

v(x)  s  (the  number  of  values  u  (x)  <  0) 

n 

Sylvester1 s  Law  implies  that 

v(x)  =  (the  number  of  J's  eigenvalues  \  < x)  • 

Therefore  any  selected  eigenvalue  can  be  computed  as  the  limit  of  a 

sequence  of  nested  intervals  [x  ,  x  ]  with 

-m  m 


and 


x  <  x  <  x  .,  <  x  for  all  m  , 
-in m+1  —  mtl  -  m  1 


x  -  x  — >  0  as  m  -» 00  . 

m  -m 


The  mechanism  by  which  the  successive  values  x  and  x  are  chosen  is 

—m  m 

of  no  consequence  here;  a  bisection  method  could  be  used  (cf,  Wilkinson 

(1962)),  though  that  is  slow.  A  faster  algorithm  has  been  produced  by  the 

author  and  Jim  Varah  (1966).  But  the  error-analysis  is  independent  of  the 

way  in  which  the  values  x  and  x  are  chosen  provided  they  have  the 

—m  m 

properties  listed  above. 

So  far  we  have  not  seen  anything  very  new.  Indeed,  the  function 
v(x)  is  just  the  number  of  variations  of  sign  in  the  Sturm  sequence 


<P0  B  >  $^(x)  )  ^2^x) 


which  has  been  in  use  for  over  a  decade  to  compute  the  eigenvalues  of 
symmetric  tri-diagonal  matrices.  (See  Wilkinson's  book  (1965)  P*  299-312. 

Also  see  Householder's  book  (1964)  p.  86-7  ex.  10  and  11,  and  p.  175  ex.  14; 
his  <pn  differs  from  ours  by  a  factor  of  (-l)n  .)  However,  the  cp-recurrence 


n-2 


takes  more  time  on  most  machines  than  does  the  u- recurrence;  and  over/ 
underflow  is  an  inescapable  complication  in  the  cp-recurrence  whereas  the 
u-recurrence  can  be  rendered  almost  immune  to  over/ underflow.  These  are 
the  reasons  Davison  gave  for  his  preference  of  the  u-recurrence.  Unfor¬ 
tunately,  he  died  before  he  had  the  chance  to  show  how  well  behaved  his 
method  could  be.  The  task  of  analysis  is  now  ours. 
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Over/ Underflow  on  the  machine 

Over/ underflow  in  the  u- recurrence  can  easily  be  rendered  insignificant 
by  a  proper  preliminary  scaling  of  the  data  and  b^  .  The  description 
of  the  scaling  process  begins  with  a  definition  of  certain  machine  constants: 

ft  is  the  greatest  floating  point  number  normally 
represented  directly  in  the  machine. 
r\  is  the  smallest  positive  (non-zero)  floating  point 
number  normally  represented  directly, 
e  is  the  smallest  positive  floating  point  number 
such  that  the  computed  value  of  1.0  +  e  differs 
from  1.0  after  it  is  rounded  or  truncated  to  the 
precision  being  carried. 

The  following  table  lists  typical  values  for  these  parameters: 


IT 


Only  a  few  machines  adhere  to  this  precision  in  all  double  precision  operations 


We  must  also  make  certain  assumptions  about  the  treatment  of  arithmetic 
over/ underflows,  because  they  will  occur.  First,  we  assume  that  whenever 
any  floating  point  arithmetic  operation  (+  #-,*,/ )  underflows  its 
result  will  be  cleared  to  zero.  Second,  we  assume  that  whenever  any  opera¬ 
tion  overflows  its  result  will  be  set  to  or  with  the  correct  sign. 

The  preservation  of  sign  after  overflow  is  essential  .•  Fortunately,  these 
conventions  for  the  treatment  of  over/ underflow  are  widely  used  on  many 
machines,  including  the  IBM  709^  and,  possibly,  the  Burroughs  B5500. 
Unfortunately,  the  new  IBM  ?60  series  hardware  forgets  the  sign  after  over¬ 
flow,  but  presumably  that  oversight  will  soon  be  corrected,  It  is  possible 
to  prevent  the  u-recurrence  from  overflowing  at  all,  but  to  do  so  costs 
a  noticeable  retardation  on  most  computers,  as  we  shall  see. 

If  over/ underflow  is  treated  as  described  above,  any  over/ underflow 

occurring  in  the  u-recurrence  will  be  practically  inconsequential  for 

reasons  to  be  given  later.  Therefore  we  must  make  a  third  assumption;  we 

assume  that  the  program  can  inhibit  the  production  of  diagnostic  over/ underflow 
messages  and  can  ignore  any  over/ 

underflow  indicators  that  might  otherwise  serve  as  superfluous  distractions 
during  the  computation  of  the  u's  .  (This  is  not  meant  to  imply  that 
those  indicators  are  superfluous  in  any  other  context.  Quite  the  contrary!) 

If  all  three  assumptions  about  the  treatment  of  over/ underflow  are 
valid  then  they  cope  with  the  problem  far  more  simply,  elegantly  and  econ¬ 
omically  than  any  other  scheme  known  to  the  author.  There  is  reason  to 
doubt  that  any  comparable  scheme  could  ever  be  devised  for  the  cp-recurrence . 
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Scaling: 


Now  let  the  scale  factor  a  be  defined  as  the  largest  power  of  the 
machine's  arithmetic  base  which  satisfies 

a | a ^ |  <,  Tfi  and  a|b^|  <  tQ  for  all  i  , 

where 

.  l/  4  - 1/  2 

The  significance  of  this  constant  t  is  that  over/ underflow  will  later  be 
shown  to  contribute  an  absolute  error  no  larger  in  magnitude  than  about 
•  lub„ ( J)  to  the  computed  eigenvalues.  Tne  values  of  t  tabulated 

O 

above  show  how  small  t  usually  is  compared  with  the  rounding  error  level 
Evidently  over/ underflow  will  hardly  ever  restrict  the  range  of  magnitudes 
spanned  by  the  accurately  computed  eigenvalues  of  J  nearly  as  much  as  do 
rounding  errors. 

Normally  a  is  approximately 


but  there  are  exceptional  cases  where  that  expression  would  overflow,  so  0 
must  be  set  instead  to  the  largest  power  of  the  machine!)  arithmetic  base. 
These  cases  are  ignored  in  what  follows  because  they  are  susceptible  to  a 
simpler  analysis  with  the  same  results  as  are  demonstrated  below. 

After  a  is  known,  the  matrix  J  is  scaled  by  being  replaced  by  oJ  . 
Since  a  is  a  power  of  the  base  there  are  no  rounding  errors.  But  under¬ 
flows  may  occur.  These  underflows  result  in  the  annihilation  of  at  most 


Tfl/tnax.  (max^  |aj  ,  max  J 


those  elements  e  and  b,  which  satisfy 

i  i  J 

laj  <  T}1/2?  lubg(j)  or  |bi|  <  ^ 2t  lubg(j)  . 

These  perturbations  are  negligible  compared  with  what  follows,  so  they  may 
be  ignored.  Later  the  computed  eigenvalues  k^  will  be  unsealed  by  dividing 
them  all  by  a  ,  Any  over/ underflow  which  occurs  here  is  fully  deserved 
and  must  be  reported  by  the *d la gnostic  machinery  mentioned  above  to  indicate 
that  some  eigenvalues  (just  the  ones  that  over/ underflow)  cannot  be  repre¬ 
sented  in  the  normal  way  without  over/ underflow,,  Nothing  more  need  be  said 
about  scalings  we  merely  assume  henceforth  that 

tH  <  lubg(j)  <  3tQ 
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Two  programs: 


Now  is  the  time  to  write  out  the  u- recurrence  explicitly  in,  say, 
FORTRAN.  There  are  two  versions,  according  as  overflow  is  prevented  or 

allowed.  Both  versions  begin  by  constructing  the  arrays  BB  and  A  con 
taining 


BB(I)  = 
A(I)  =  a-j. 


for  I  =  1,2,  . .  N  . 


If  <  V^T  then  BB(l)  will  underflow  to  zero,  but  this  amounts  to 

a  perturbation  of  no  more  than 


yr  -  t  •  (ta)  <  t  •  iubg(j) 

in  the  given  matrix  J  ,  and  is  included  in  the  error  analysis  given  later. 
Note  that  BB(l)  =  0.0  by  definition,  and  that  we  can  assume  that 

|x|  <  lubE(j)  < 

is  satisfied  by  any  number  X  which  might  usefully  be  considered  as  an 
estimate  of  an  eigenvalue. 

Here  is  the  segment  of  code  which  prevents  any  overflow  in  the  u- recur¬ 
rence  ;  the  constant  RTETA  has  the  value 


RTETA  -  yr  • 
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The  FORTRAN  symbol  ".GT." 


stands  for  ">"  . 


U  =  1.0 

NU  =  0 

DO  3  I  =  1,N' 

TJ  *  (A(I)  -  BB(I)/U)  -  X 
IF  {U  .GT.  RTETA)  GO  TO  3 

1  IF  (U  .GT.  -  RTETA)  U  =  -  RTETA 

2  NU  =  NU  +  1  ...  when  U  <  0  . 

3  CONTINUE 

. ..  Now  NU  =  v(x)  . .. 

,  it 

bi/ui 


Whenever  the  computed  value  of  u^  lies  between  -y?  and 
is  replaced  in  statement  1  by  -’yfi]’ .  Consequently  the  quotient 
never  exceeds 


(Tfi)2/yy  =  n  , 


so  overflow  is  impossible.  Of  course,  u^  may  have  been  decreased  in 
statement  1  by  as  much  as  ,  but  this  too  is  no  larger  than  might 

have  been  caused  by  decreasing  by  the  allowable  perturbation 

2~\p\  <  2t  •  lubg(J) 

Here  is  a  simpler  and  faster  program  segment  which  is  useable  whenever 
overflow  is  treated  according  to  the  conventions  described  above.  The 
constant  ETA  -  r\ 
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U  =  1.0 


NU  =  0 


DO  3  I  =  1,N 

U  =  (A(I)  -  BB(l )/U)  -  X 


IF  (U)  2,  1,  J 

1 

U  =  -  ETA 

if 

U  was 

0  . 

2 

NU  =  NU  +  1 

if 

U  was 

<  0  . 

3 

CONTINUE 

...  Now  NU  =  v (X)  , . . 

Note  that  whenever  any  vanishes  it  is  replaced  in  statement  1  by 
u^j.  «  -tj  to  forestall  a  subsequent  division  by  zero.  Whenever  (rarely) 
any  u^.  overflows,  its  sign  remains  unchanged  so  that  NU  is  treated 
correctly;  then  uI+1  is  in  error  because  the  computed  value  of  bj/uj 
must  be  larger  in  magnitude  than  it  should  be.  But  the  error  in  uI+1 
is  no  worse  than  might  have  been  caused  by  perturbing  a^+^  by  at  worst 


(tn)2/n  =  t(tA)  <  t  •  lubq(j) 


s 


And  underflow,  if  it  occurs,  causes  no  more  perturbation  than  r|  ,  which 
is  negligible. 

All  told,  these  subterfuges  for  circumventing  the  ill  effects  of  over/ 
underflow  cause  the  computed  value  NU  to  be,  instead  of  v(X)  ,  some 
value  that  would  have  been  obtained  had  J  first  been  changed  in  each 
element  by  at  most 

2t  *  lubG(j)  (in  the  first  program)  or 


2k- 


t  0  lubg(j)  (in  the  second  program) 

before  v(X)  was  computed  without  any  intervention  on  behalf  of  over/ under¬ 
flow.  These  perturbations  will  be  shown  later  to  affect  the  computed 
eigenvalues  by  no  more  than  4t  «  lubg(j)  .  First  we  should  consider  one 
last  programming  detail. 

The  k  eigenvalue  is  the  k  jump-point  of  the  integer  valued 

function  v(x)  ; 


lim  v(x)  <  k  <  lim  v(x) 

X-^k,  “  X-*d.  + 

k  k 

(The  multiplicity  of  the  jump-point  is  just  the  difference  between 

the  upper  and  lower  limits.)  Can  a  similar  statement  be  made  about  the 
computed  approximation  NU(x)  ?  If  so,  any  algorithm  that  works  properly 
for  the  exact  function  v(x)  will  work  properly  for  its  approximation 
NU(X)  .  If  not,  if  NU(x)  could  have  more  than  N  jump-points,  then 
great  care  would  be  required  to  design  the  algorithm  in  such  a  way  that 
it  could  not  be  confused  by  spurious  jumps  down.  As  it  happens,  no  such 
care  is  required  on  most  machines. 

We  shall  demonstrate  below  that,  despite  rounding  errors  and  over/ 
underflow,  the  computed  function  NU(x)  is  a  monotonic  non-decreasing 
integer  valued  function  of  X  with  just  N  jumps.  The  only  assumption 
is  that  each  arithmetic  operation  executed  by  the  machine  is  a  monotonic 
function  of  its  arguments  despite  rounding.  For  example,  if  A  ,  B  and 
C  are  all  positive  numbers  represented  in  the  machine,  and  if  the  FORTRAN 
program 
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XI  = 


A  +  B 


Y1  =  (A  +  C)  +  B 
X2  =  A  -  B 
Y2  =  (A  +  C)  -  B 
X3  =  A  *  B 
Y5  =  (A  +  C)  *  B 
X4  =  A  /  B 
Yk-  =  (A  +  C)  /  B 
X5  --  B  /  (A  +  C) 

Y5  =  B  /  A 

is  executed,  then  XI  <  YI  for  all  I  =  1,2,3#^  or  5  .  This  assumption 
is  certainly  valid  for  single  precision  computations  on  all  of  the  machines 
listed  in  the  table  above.  Indeed,  the  builder  of  any  machine  which  failed 
to  satisfy  this  assumption  should  be  ashamed  of  himself. 
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The  Monotonicity  of  NU(X): 

The  monotonicity  of  NU(x)  will  be  derived  as  8  consequence  of  the 
properties  of  the  successive  values  of  U  ,  for  which  some  notation  is 
required.  Given  any  argument  X  ,  a  number  representable  in  the  machine, 
either  program  above  will  produce  a  sequence  of  values  U  (x)  and  NU^(X ) 
the  values  taken  by  U  and  NU  respectively  after  statement  5  has  been 
executed  for  the  n"^  time.  In  particular, 

U  (X.)  *  1.0  and  NU  (X)  =  0  ; 

o  o 

U-^(X)  --  L A ( 1 )  -  X]  rounded  etc.  and 
NU,(X)  =0  if  X  <  A(l)  -  6 

=  1  .if  X  >  A(l)  -  6  , 

where  0  ■v.'Tf  in  the  first  program 

=-•  r\  in  the  second  program  . 

Note  that  no  U  (X)  can  lie  closer  to  zero  then  +0  or  -0  .  At  the  end 
n 

of  the  D@-loop, 

U  *  UN(X)  and  NU  *  NU(X)  -  NUN(X)  . 

The  interesting  values  of  X  arc  those  where  some  Un(x)  changes 
sign.  These  points  shall  be  identified  precisely  with  the  aid  of  a  not at. it 
X/  for  the  successor  of  X  ;  if  X  is  a  number  representable  in  the  ma¬ 
chine  and  eligible  to  be  an  argument  for  the  programs  above,  then  X/  is 
the  next  larger  eligible  argument.  Normally  X'  will  exceed  X  by  one 


unit  in  their  last  place  being  carried  in  the  computation. 

A  "zero"  Z  of  Un (X )  is  now  defined  to  be  any  argument  Z  which 
satisfies  both 


U  (Z)  >  6  and  -e  >  U  (z')  . 

n  -  —  n' 

A  "pole"  Y  of  U^(X)  is  any  argument  Y  which  satisfies 

U  (Y)  <  U  (Y;)  . 

nx  '  ns  ' 

Un(X)  can  change  sign  only  at  a  zero  or  a  pole,  though  U  (X)  may  fail 
to  change  sign  at  some  poles.  Between  any  two  zeros  of  U  (X)  must  lie 
at  least  one  pole  where  Un(X)  changes  sign,  and  possibly  some  other 
poles  where  Un(x)  does  not  change  sign.  Let  us  examine  these  poles  more 
closely. 

If  Y  is  a  pole  of  Un(X)  then 

rBB(n)/Un-1(Y)]  >  [BB(n)/Un->1(Y/)] 

because  the  contrary  relation  would  prevent 

U  (X)  =  [[A(n)  -  [BB(n)/U  _  (X)]  ]  -  X]  , 
n  n-l 

where  each  pair  of  brackets  means 

"round  [...]  and  take  care  of  over/ underflow,  if  any"  , 
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from  increasing  when  X  moves  from  Y  to  Y/  .  (Note  that  the  over/ 
underflow  subterfuges  do  not  destroy  the  monotonicity  of  the  arithmetic 
operations  even  if  U  (x)  has  to  be  replaced  by  -•©■.)  Therefore  either 
Y  is  a  pole  of  ^(X)  where  Un  ^  does  not  change  sign,  or  Y  is  a 
aero  of  ^(X)  .  A  backward  induction  yields  the  following  statement: 

If  Y  is  a  pole  of  U^(X)  ,  then  there  exists 
some  positive  integer  m  <  n  such  that 

U  (Y)  >  0  >  U  (Y')  , 

mN  mx  '  • 

and  for  all  integers  i  (if  any)  strictly  between 
m  and  n 

UX(Y)  <  U^l') 

with  no  change  of  sign. 

We  ab serve  that  U^(x)  has  no  poles  and  just  one  zero.  Therefore, 

as  U_(X)  is  carried  from  U0(-ft)  -  ft  to  IL  (ft)  =  -ft  ,  it  can  have  at 
d  d  d 

most  two  zeros  separated  by  one  pole  where  changes  sign,  or  one  zero 

and  one  pole  where  does  not  change  sign,  or  one  zero  and  no  poles  if 

d 

BB(2)  is  very  tiny.  In  all  cases  one  can  verify  with  ease  that  NIL ( X ) 

d 

is  a  monotonic  non-decreasing  function  of  X  with  at  most  two  distinct 

jumps  from  NIL  (-.ft)  -  0  to  NIL  (ft)  -  2  .  Rather  than  extend  this  desired 

d  d 

property  to  NUn(x)  for  all  n  by  a  long  constructive  argument,  we  shall 

show  that  a  failure  of  NU  (X)  to  be  monotonic  would  create  a  contradiction. 

n 

ft  9 


Let  n  be  the  smallest  integer  for  which  NU  (X)  is  not  a  monotonic 
function.  Obviously  n  >  1  .  Suppose  NUn(X)  falls  to  be  monotonic  at  Y 
since 


0  «  NU  (-n)  <  NU  (X)  <  n  o  NU  (ft)  , 

n  —  n  ~  n 

the  failure  must  take  the  form 

KUn(y)  >  NUn(Y')  . 

However*  our  hypothesis  about  n  imp3.ies  that  NU  ^(x)  is  monotonic* 
which  means 


Also, 

NU  (X)  -  NU  ,  (X)  -  0  if  U  (X)  >  0 
n  n-i.  n 

»  1  if  Un(X)  <  0  , 

so 

0  >  NUn(Y')  -  NU n(Y) 

"  {nub(y')  -  NUn_J ( y' )) 

+  •  NUn-l(Y» 

+  {NU  ,(Y)  -  NU  (Y)) 
n-J.  n 

>  Co]  +  {0}  f  (-1]  «  -l  . 
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But  this  Implies,  [term]  by  {term}  ,  that 

NU  (Y')  -  NU  .  (y')  and  U  (y' )  >  0  , 
n  n-1  n  * 

NU  .  (y' )  -  NU  _  (Y)  «=  0  ,  and 

n-1  n-1'  ' 

NU  (Y)  «  NU  , (Y)  +  1  and  U  (Y)  <  0  . 

n'  '  n-1'  1  n'  1 

Evidently  Y  is  a  pole  of  U^(X )  .  Therefore  there  exists  some  positive 

integer  m  <  n  for  which  Y  is  a  zero  of  Um(X.)  }  wo  shall  have 

U  (Y)  >  0  >  U  (y' )  . 
m'  m'  ' 

Therefore  NU  (Y;)  «  NU  .(y')  +  1  because  U  (y')  <  0  , 

m  '  m-1  m  ; 

>  NU  _2.(Y)  +  1  by  monotonicity  , 

m  NU  (Y)  +  1  because  U  (Y)  >  0  . 

m  m 

Also,  if  there  arc  tiny  integers  :1  strictly  between  m  and  n-1  , 

NU^Y')  -  NU^tt')  ■  NUi(Y)  -  NU1ib1(y) 
because  Y  is  u  pole  of  U^  with  no  change  in  sign,  Therefore 

nu  ,(y')  -  nu  Ay)  -  nu(y')  -  nu  (y)  >  i.  , 

n-1'  .  '  n-1  m'  m  — 

who rou s  wo  saw  above  that 


NU  (Y')  -  NU  .  Of)  »  0 
n  -  .1,  n  -  .i 


d.l 


This  contradiction  proves  that  NU(X)  is  a  monotonic  non-decreasing  function 
of  X  ,  as  desired. 

It  seems  surprising  that  so  strong  a  result  can  be  proved  with  no 
appeal  to  the  continuum,  nor  any  estimate  for  the  errors  in  the  values  Un(x) 
On  the  contrary,  the  values  of  Un(X)  can  be  completely  different  from  the 
mathematically  exact  values  un  that  would  have  been  obtained  without 
rounding  errors  nor  over/ underflow,  even  to  the  extent  of  having  the  wrong 
signs.  Fortunately,  the  errors  in  the  intermediate  results  U  (X)  are 
of  no  interest  beyond  an  assurance  that  the  errors  are  not  haphazard.  And 
the  behaviour  of  NU(X)  provides  just  that  assurance. 


Bounding  Rounding;  Errors; 


The  next  step  is  to  show  that  if  |x|  <  iub  (j)  then  NU(X)  is  pre- 

lij 

cisely  the  value  that  v(X)  would  hove  taken  if  J  had  been  replaced  by 
some  nearby  matrix  J(X)  and  all  computations  hod  been  carried  out  in¬ 
finitely  precisely  with  neither  rounding  errors  nor  over/ underflow  subter¬ 
fuges.  The  principles  behind  the  analysis  that  follows  ore  very  much  like 
those  to  be  found  in  Wilkinson's  books  (I96j5,  19 65).  We  shall  try  to  describe 
the  elements  of  J(X')  in  terms  of  the  numbers  that  actually  appear  in  the 
arithmetic  registers  of  the  machine  during  the  computation,  and  in  terms  of 
the  rounding  error  bound  e  tabulated  above  for  several  machines.  The 
ideas  involved  are  best  illustrated  by  the  following  examples. 

The  FORTRAN  assignment  statement 

C  n  A  *  13 

will  not  replace  C  by  the  product  of  A  and  B  }  but  will  instead  sot 
C  to  u  value 


C  *  (,l  +  7  )AB 

in  which  | 7 |  is  normally  bounded  by,  soy, 


<  Ci 


t 


Note  that  A  ,  B  and  the  now  value  C  ure  defined  quite  precisely,  and 
satisfy  the  previous  equation  exactly.  The  only  unknown  quantity  is  7  , 


but  | 7 |  Is  bounded  by  a  known  value  e  except  when  over/ underflow  inter¬ 
venes.  Similarly,  the  assignment  statement 


actually  stores  a  value 


where 


c  =  a/b 


C  «  (1  +  3)A/B 

|p|  <  c  . 


As  a  matter  of  fact,  the  situation  Is  not  always  as  described  above. 
In  double-precision  the  values  of  3  and  y  can  be  as  large  as  on  a 

7094,  on  a  B5500,  and  l6e  on  a  360,  These  unnecessarily  large 

errors  are  so  repugnant  to  the  author  that  he  takes  the  .liberty  of  passing 
them  directly  from  the  machines'  manufacturers  to  the  reader,  who  may  ac¬ 
commodate  them  by  multiplying  a  in  the  bounds  given  below  by  whichever 
factor  is  appropriate  for  his  machine.  For  a  similar  reason,  the  author 
chooses  to  presume  that  the  FORTRAN  statement 


C  *A  +  B 


causes  C  to  take  precisely  the  value 

C  «  (A  +  B)/(,l  +  a) 

with  | ex |  <  c  .  Only  in  double -precision,  end.  then  only  on  some  machines, 
is  it  necessary  to  replace  the  last  two  relations  by 

C  «  (1  +  a)A  +  (1  1  3)B 
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but  this  will  weaken  the  error  bounds  to 


with  | o; |  <  €  and  |p|  <  e  , 
be  given  below  by  a  factor  no  larger  than  two. 


The  Construction  of  J(X)  near  J: 

The  first  step  in  the  construction  of  J(x)  is  the  definition  of 
certain  values  and  corresponding  respectively  to  a^  and  b^  . 

Let  us  set 


p 

B^  »  bj  if  Bj  does  not  underflow  , 
a  0  if  Ibjl  <  yT  • 

In  either  event,  |b^  -  b^|  <  ^rf  .  And 


BB  B  BB  ( I )  «  (1+  , 

where  ^  is  the  relative  error  due  to  multiplication  and  satisfies 
If,;*  |  <  g  .  (The  primes  used  during  the  construction  of  J(X)  do  not  de- 
note  successors.) 

The  value  of  A^  depends  upon  X  ,  arid  differs  from  only  to  the 
oxtont  required  to  compensate  for  the  effects  of  over/ underflow.  Res  sons 


have  already  been  given  why  we  should  expect  that 


K  -  ojl  <  2  Vn”  in  the  first  program  , 

<y^  in  the  second  program  , 

execpt  for  an  ignorable  contribution  no  larger  than  r\  .  Let  us  stay  with 
the  second  program  from  now  on,  and  ignore  not  only  the  scaling  underflow 
error  r\  by  setting 


A(I)  5  Sj  , 

but  also  agree  to  ignore  the  comparable  error  induced  by  underflow  or  state¬ 
ment  1  during  the  computation  of  Uj  . 

The  dissection  of  the  FORTRAN  statement 

U  «  ( (A(l)  -  DB(I)/U)  -  X) 

to  find  its  rounding  errors  is  an  inductive  process.  For  1  »  1  we  define 

A^  w  A(l)  *  und 

v.|  »  A.j  -  X  precisely,  nnd 

U1  -  f  [ A ( 1. )  -  0.  ]  -  X]  «  Vj/d  +  o£)  , 


where 


|a!jj  <  c.  ,  Evidently  sign(v^)  olgn(U.j ) 


—  1  . 


The  induction  hypothesis  is  that 


We  also  set  a'  =  p''  =  0  and  v 
11  o 

for  n  =  1,2,...,  1-1  we  may  write 


v  =  A  -  (1  +  p" )(l  +  a  n  )(l 
n  n  '  n'  v  n-1' ' 


a'  )BB  / v  n  - 
n-L  n  n-1 


(1  +  a')X  , 
n  ’ 


lAn  -  \\  <  VT  - 


(pj  <  €  '  lanl  <  €  '  l°vJ  <  €  ' 


and 


either  v  «  (l  +  0?  )(l  +  a1  )U 
n  n  '  n'  n 


or 


lu  I  «  n  and  0  <  U  /v  <1  and  a'  »  d'  <=  0  . 

1  n1  n'  n  n  n 


In  any  case,  aign(vn)  »  sign(U^)  />  0  . 


The  hypothesis  is  obviously  true  for  n  «  1  ,  since  BB^  «  0  .  Note  that 

U  represents  what  wen  enplier  referred  to  as  U  (X)  ,  and  is  a  number 

actually  stored  in  the  computer,  The  values  BB  and  a  ore  also  stored 

n  n 

in  the  computer,  but  will  not  be  stored  if  it  differs  from  a^  ,  nnd 

v  in  u  figment  of  the  imagination  except,  for' its  sign, 
n 

Now  for  the  advance  to  n  *  1  ,  The  first  value  to  be  inspected  is 


[  BB ( 1 )/ Uj _ ] rounded  »  { 1  +  p."  )^1/^lml  , 


where  | (VJ, |  <  t  unless  overflow  occurs.  (Underflow  is  being  ignored.) 

If  this  quotient  overflows  then  the  remaining  arithmetic  operations  are 
Irrelevant  because  the  scaling  has  ensured  that  neither  A(l)  nor  X  can 
be  bigger  thnn 


3tfi  «  efi  ; 


therefore  overflow  will  cause  to  be  given  the  value 

Uj  =  -n  sign(UI-;i)  , 


and  we  may  define 


Aj  s  ar  ,  p*  3  a'  b  a"  3  o  ,  and 
VI  3  <A!  -  “A-l*  *  X  • 

These  values  satisfy  the  induction  hypothesis. 

If  the  quotient  [BB^/U^]  does  not  overflow,  there  is  still  the 
possibility  that  the  previous  quotient  £0Bi..i/UI-2^  overflowed  to  be 
considered.  In  this  case  |uX-1J  B  ft  and  0  <  <  1  ,  and  we 

define 

Aj  *  A(I)  -  [I3B(I)/UI-;l]  +  (1  +  Pj)BB1/v1.-1  . 

Evidently 

|AI  -  ox|  -  (1  +  p"  )BBI|l/UI_1  -  1/Vj.j.I 

<  (;l  +  «)(tn)s|i/n|  »*  (i  +  <■)  VT  • 


The  factor  1  e  is  unimportant  and  shall  be  dropped.  Note  too  that 


A(I)  -  [BB(I)/UI-;l]  «  Ax  -  (1  +  Px)BBI/vI-;1 

'  Ai  ■  (l  +  p'xH1  +  +  ai-.i)BVvi-i 
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The  la at  relation  is  satisfied  too  if  neither  [BB(l)/Uj  nor 
[BB(T.-1)/U^  n]  overflowed,  in  which  case  we  set  A  .  s  A(l)  -  ,  because 

then  vj  “  (l  +  c/j  ^  )  (l  +  Qfj  ^  . 

Continuing,  if  [BB(l)/lI^  ^ ]  does  not  overflow  then  the  value  stored 
for  Uj  will  be 


U-j.  **  [[At  -  (1  +  f\  ) (l  +  +  a'_1)BBI/v];_I]  -  X] 

«  ((Aj  -  (1  +  pj)(l  f  o^Kl  +  a'-1)BBI/vI_1)/(l  +  a') 

-  X)/(l  +  c/') 
s  vI/((l  +  a")(l  +  a' ) ) 


where  the  rounding  errors  of  addition  are  bounded  by  |a£|  <  e  and 
|c/'. |  C  e  .  Thio  result  advances  the  induction  from  n  *  I  -  1  to  n  »  I 
os  do aired,  and  lays  a  firm  foundation  for  on  error  bound  for  the  eigen- 
va lues . 


Lot  the  matrix  J(X)  be  defined  now  to  hove 


A  -  a;X  in  place  of  o  ,  and 
n  n  1  n  ; 


'")(!  +  </  ,)(l  +  a'  )  B  . 
n/N  n-1  n-1  n-1 


in  place  of  b  , 
n  *  i 


Certainly  J(X)  is  clone  to  J  ;  more  precisely,  but  neglecting  terms 

p 

of  order  <:  and  i|  , 


|  J (X)  -  j|  <  Pc-:  |  J  -  diag  j|  »•  a|x|l  +^/x\  H 


e  lemon  tw  .1.  so ,  who  re 


J  -  diag  J  3  only  the  off-diagonal  terms  b^  in  J 
H  -  the  tri-diagonal  matrix  with  all  elements  =  1 
Denote  the  eigenvalues  of  J(X)  by 

<  V>(X)  <  ...  <  ^(X) 

to  correspond  with  the  eigenvalues  k^  of  J  , 


ko 


The  Absolute  Error  in  K.  : 

K 


The  reason  for  constructing  J(x)  was  that  NU(x)  would  be  the  number 
of  J(x)'s  eigenvalues  \^(X)  <  X  ,  and  now  this  can  be  proved.  For 


NU(X)  =  the  number  of  values  U  (X)  <  0 

«  the  number  of  values  v  <  0  , 

n  ’ 


and  the  v  are  to  J(x)  what  the  uncontaminated  values  u  are  to  J  . 
n  n 

And  each  eigenvalue  \^(X)  differs  from  the  corresponding  \  by  no  more 


than 


lub„(j(X)  -  J)  <  2e  lub,,(|j  -  diag  j|  )  +  e|x| 

u  “  o 

+yr  iub0(H)  . 


Here 


V?  lubg(H)  <  jVT  ■  lubs(J)  , 


2e  lub„( | J  -  diag  j| )  <  Be  lub0(j) 


by  virtue  of  tho  more  general  form  of  our  earlier  result  no.  2  with 


la2  +  fi2  *“\A)  t  (2c;  )2 


Finally,  the  only  values  of  X  that  will 


concern  us  below  are  those  which  approximate  some  eigenvalue  ,  so 
we  can  certainly  assume  that  e|x|  <  e  max,^ | |  »  c  lub0(j)  to  within 


a  negligible  extra  error  of  the  order  of  c  |X|  .  For  those  values  of  X 


MX)  -  <  r  ■  5(e  +  r)max  |\ 

J  |J  u 


we  have 


as  a  bound  for  the  difference  between  the  corresponding  eigenvalues  of 
j(X)  and  of  J  .  This  means  that,  as  X  varies  over  the  allowable 
arguments,  each  eigenvalue  \^(X)  remains  confined  to  some  fixed  interval 

\i  -  r  <  \t(X)  <  \i  +  r  , 

We  have  already  seen  that  for  any  given  k  there  is  precisely  one  value  X 
which,  with  its  successor  X^  ,  satisfies 

NU(Xk)  <  k  <  NU(Xk)  ) 

these  values  can  easily  be  computed,  And  the  relationship  between  NU(x) 
and  the  \^(X)  telle  us  that 


Xk  -  W  -  \  +  r  and 


Since 

0  <  X^  -  Xk<  2e  maxJ\J  on  u  rounding  muohino 
<  c  on  a  truncating  machine 

(we  might  os  well  uimuma  now  that  arithmetic!  :!,«  rounded),  we  con  accept 
either  'Xk  or  Xk  as  un  approximation  to  \k  and  commit  an  error  no 
larger  than 


Thin  bound  compares  favour ably  with  that  obtained  by  Wilkinson  ( 1965 > 
p.  *,04-5)  tor  the  Sturm- sequence  (cp-recurrence )  algorithm  in  the  absence 
of  over/ underflow.  In  our  notation  his  bound  is  17 e  maXj)\^|  ,  although 
the  use  of  our  more  refined  methods  reduces  this  to  8.75€  maXj|\j|  . 

This  bound  is  not  appreciably  increased  if  Wilkinson's  1962  program  is 
amended  to  cope  with  over/ underflow,,  but  then  the  cp» recurrence  becomes 
much  slower  than  the  u-recurrence .  Therefore  the  u-recurrence  has  all  the 
advantages  of  speed,  simplicity  and  accuracy  over  the  tp-reeurrence .  On  a 
computer  using  rounded  binary  floating  point  arithmetic,  the  biggest  eigen¬ 
value  c/an  be  computed  to  within  a  guaranteed  relative  error  of  2.5  units 
in  its  last  place,  and  no  eigenvalue  will  suffer  a  larger  absolute  error. 
For  chopped  arithmetic  the  guarantee  is  4  units  in  the  lost  place.  These 
bounds  are  impressively  small;  but  they  are  substantially  larger  than  most 
of  the  errors  observed  in  practice. 

Why? 


4'i 


The  Clusa  of  neighbours  of  J  t 


A  nicer  appreciation  of  the  accuracy  of  the  u-reeurrence  can  be  achieved 
through  the  consideration  of  the  class  ^  of  symmetric  tri-diagonal  matrices 
v5r  which  satisfy 


|j  «  <  2e|j  -  ding  j|  . 

For  example,  on  the  rounding  binary  computer  mentioned  above  the  class  ^ 
consists  of  those  matrices  J  obtained  from  -T  by  changing  each  off- 
diagonal  element  of  J  by  at  moat  one  unit  in  its  last  place,  The  set 


9 


is  a  convex  set  In  the  sense  that  if  j  and  J.  are  members  of 

u  1 


$ 


then  so  are  all  matrices  of  the  form 


tJx  4*  (1  -  t)JQ  for  0  <  t  <  1 


lying  "between1'  ifQ  an l  Jj  ,  Much  matrix  t  in  ^  has  a  set  of  eigen- 
valued 


1  “  d  —  “  N-i  —  N 


and  as  J  varies  over  ^  each  eigenvalue  ^  varies  over  some  set 
wnich  can  also  bo  shown  to  be  a  closed  convex  set,  In  other  words,  ouuoata- 
ted  with  the  class  ^  of  matrices  j  in  the  set  of  N  intervals 

^  the  cot  of  nil  x  «  for  some  j  in  ^ 
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cun 


Home  o.f  these  intervale,  may  overlap,  but  it  is  soon  seen  that  no  A 

K 

be  contained  strictly  inn id e  another.  Therefore  the  intervals  A,  shore 

k 

the  same  ordering  an  the  eigenvalues  ..  Obviously  Aj  is  contained  in 


"k 


ile  rr.'iXj  \  k^  j  <  x  <  +  2e  imlXj|\j| 


but  the  interval  A^  hardly  ever  oooupi««  more  than  a  small  fraction  of 
that  latter  interval. 


The  significance  of  the  interval  A^  is  that  for  most  practical  pur¬ 
poses  any  number  in  A^  la  as  acceptable  an  approximation  to  as  any 
other.  Ouch  might  be  the  case,  for  example,  if  each  off-diagonal  element 
b^  of  J  w^re  independently  in  error  by  as  much  as  5>e|b,  |  because  of 
previous  rounding  errors.  The  independence  of  the  errors  is  essential; 
correlations  among  the  errors  in  the  bj  could  conceivably  cause  the 
eigenvalues  of  0  not  to  be  in  error  at  all,  as  would  be  the  cane  if 


were  erroneously  computed  us 


As  long  uu  the  errors  in  the  bi  are  independent,  the  width  of  the  inter¬ 
val  A^  is  un  indication  of  the  extent  to  which  must  bo  regarded  aw 

inoscs pnbly  uncertain.  And  that  the  error  introduced  by  the  programs 
analyzed  hero  contributes  negligibly  to  this  inescapable  uncertainty  uhull 
now  be  demonstrated. 


45 


Let  J(X)  be  defined  to  have  elements  4  *  a  and 

N  n  n 

^n-i  +  f^Kl  +  P*)(l  +  O^Ml  +  o^)  bn-]  ,  where  the  Greek 

letters  were  defined  during  the  construction  of  J(X)  .  Except  for  terms 
2 

of  order  ts  which  shall  bo  ignored, 


-  b  1 |  <  |b  1 | 

1  n-1  n-1 '  -  1  n-1' 

so  J(X)  belongs  to  J)  and  each  eigenvalue  £^(x)  of  J(X)  lies  in  ito 
corresponding  interval  ,  Also, 

|i(X)  -  J(X)|  <  e|x|l  H 

except  for  ignorable  terms  of  order  r\  ,  so 

|yx)  -  \k(X)|  <  r(X)  ■  a  |  X  |  +  3'r  max,jJ\J 

And  if  X^  and  ita  uuoowwior  are  defined,  au  before,  by 

»U(Xk)  <  k  <  NU(x') 

then 

Xk  <  \(Xk)  <  <k(Xk)  +  r(Xk)  und 

W  ■  *<*£>  s s  K  ■ 

If  arithmetic  in  rounded,  then 

Xk  -  Xk<  2t  MX.(|xkl  ,  |X'|) 
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(except  possibly  when  rj  >  which  is  ignored),  Putting  ell  the 

facts  together  shows  that  either  X^  or  X*  differs  from  >t(X^)  or 
rt(Xk)  respectively  by  no  more  then 


i'c  max .  ( I  XR|  ,  |xk|)  +  3  i  naxj  |\,  |  , 

On  s  binary  machine  with  rounded  floating  point  arithmetic  we  may  summarise 
this  result  os  follows} 

The  computed  approximation  to  the  k  eigenvalue  of  J  need 

not  differ  by  more  than  1  -  units  in  its  last  place  from  the  k^1  eigen- 
value  of  some  matrix  ^  each  element  of  which  differs  from  the  cor¬ 
responding  element  of  J  by  at  most  one  unit  in  its  last  place,  plua  an 
absolute  error  of 


3*r  max.l|v,| 


from  over/ underflow  subterfuges, 

In  other  words,  if  J  is  already  uncertain  in  each  element  by  several 
units  in  ltd  lust  place,  and  if 


Ki/^iKl > 

“10 

(which  it)  not  often  a  restriction  since  jVi/e  <  10  on  each  of  the  machines 
tabulated,  even  in  double-precision),  then  the  additional  uncertainty  intro¬ 
duced  by  the  computation  of  will  be  insignificant  when  compared  with 
the  Intrinsic  uncertainty  in  caused  by  uncertainty  in  J  .  If 
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ie  intrinaioally  uncertain  by  only  a  few  units  in  its  last  place,  then  the 
approximation  to  will  be  accurate  to  within  a  few  units  in  its  last 
place  too  deapite  the  fact  that  is  much  smaller  than  max.jxj  .  This 
partially  explains  why  some  of  the  very  small  eigenvalues  of  symmetric  tri¬ 
diagonal  matrices  have  been  computed  to  such  unexpectedly  fine  relative 
precision  by  the  u-reourrence. 
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Insensitive  Eigenvalues  i 


But  why  a re  the  smaller  eigenvalues  of  J  so  frequently  (but  not 
always )  so  much  less  sensitive  to  small  relative  perturbations  in  J  than 
might  be  suggested  by  simple  examples  like  the  following? 


J  « 


Ky  «  e  ,  »  2-e 


Unlike  this  example  are  many  others  where  even  the  tiniest  eigenvalues  suffer 
relative  (rather  than  absolute)  displacements  which  are  of  a  comparable  order 
of  magnitude  with  the  relative  changes  in  the  off-diagonal  elements  of  J  , 
(Although  not  always  easy  to  explain,  it  is  often  observed  that  J's  eigen¬ 
values  are  less  sensitive  to  relative  perturbations  in  the  off-diagonal 
elements  than  to  comparable  relative  perturbations  in  the  diagonal  elements.) 
An  extreme  example  of  this  phenomenon  is  provided,  by  those  matrices  J 
whose  diagonal  elements  oil  vanish.  These  matrices  turn  up  during  certain 
computet  ion  a  of  singular  values,;  see  Golub  and  Kalian  (1965)  p.  213.  The 
methods  used  above  can  be  exploited  to  prove  that 


If  J  1b  an  NXN  symmetric  t ri**d logons  1  matrix, 
if  diug  J  ■  0  ,  and  if  |ftj|  <  t: | J |  , 


then  the  ordered  eigenvalues  of  J  and 

of  J  -1-  BJ  satisfy 

\?>\\  <  Nf.|\i|/(l  -  Ne)  , 


provided  Ne  <  1  . 
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An  outline  of  the  proof  follows.  Write 


hi  +  Bb.^  *  b^l  +  p^)  with  IpJ  <  e  . 

Without  loss  of  generality  we  may  assume 

bi  /  0  for  1  <  i  <  N  , 

Corresponding  to  the  u«reourrence  applied  to  J  ••  xl  is  the  corresponding 

v-reourranee,  say,  that  belongs  to  (J  +  BJ)  ••  xl  }  they  can  best  be  compered 

when  written  side  by  side  thus? 

U1  *  -*  i 

ui  “  *x  “  bi-?/ui«l  * 

It  is  well  known  and  easy  to  show  that  both  J  and  J  +  BJ  have  only  sim¬ 
ple  eigenvalues  and  respectively,  and  that  m  0  if  and 

only  if  x  is  an  eigenvalue  of  J  ,  and  that  vN  -  0  if  and  only  if  x 
la  an  eigenvalue  of  J  +  BJ  ,  (of,  Wilkinson  (1965)  p.  JOO.) 

4*  Vi 

Our  object  now  is  to  show  that  each  is  the  k  eigenvalue  of 

some  matrix  which  differs  from  J  +  bj  by  terms  of  order  e:X,^  rather 
than  .  There  ere  two  cases  according  as  N  is  odd  or  even. 

If  N  w  2n-l  we  define  the  factors  (l  +  y^)  via 

1  +  7n  H  1 

1  +  7iU  M  < 1  +  fV'Vd  +  7t)  ^  1  >  n 

1  +  'l-i  *  <’■  +  Pi.i^/O-  +  7*)  H'  i<n  . 


Vj 


^  +  pi-l^bi-l/vi-l 
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and  ob nerve  that 


wi  s  0-  +  ^)u|  ond  xj  a  * 


wl.  “  xl 


-  x 


wi.  "  xi 


-  X  -  (1  + 


Thin  is  juat  the  w-reourrenoe,  say,  belonging  to  the  matrix 


J(x)  -  xl  ■  J  +  bJ  +  dieg(x^)  -  xl  , 

f'iince  |J  +  ftJ  -  J(x)|  «  |x|diag|>'J 

<  | x |  dieg((l  *  0“gn+S  -  I) 

<  N«|x|l/(1  -  Ne) 

if  Ne  <  1  ,  the  kth  eigenvalue  \k  +  of  J1  *  bJ  differs  from  the 
eigonvulue  of  J(x)  by  no  more  then  Ne|x|/(3,  -  N«)  .  Put  the 
eigenvalue  of  J(\R)  1b  Just  since  oign(wj)  ■  elgn(ut)  for  all  x 
and  wN  «  Ujj  «*  0  for  x  -  .  Therefore 

Uk  +  -  \\  <  Nc|\lll/(1  -  N«) 


a  a  claimed. 
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m  sinuj-ar  scnema  worms  wnen  N  is  even.  Thus,  one  can  hardly  he 
surprised  in  this  case  when  each  computed  eigenvalue  of  such  a  matrix  J 
is  correct  to  within  N  units  in  its  last  place  despite  a  wide  variation 
in  the  orders  of  magnitudes  of  the  eigenvalues. 

The  possible  persistence  of  high  relative  precision  in  many  of  the 
tiny  eigenvalues  of  wider  classes  of  matrices  J  awaits  a  systematic 
explanation  with  predictive  powers,  in  the  absence  of  which  it  is  hard  to 
sty  when  a  small  computed  eigenvalue  has  higher  relative  precision  than  is 
implied  by  the  absolute  error  bound 

(9«  +  3t)maxj|\^|  . 


Conolutiloni 

There  are  faster  programs  than  those  described  hero,  but  none  more 
elegant  nor  more  accurate. 
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Having  estubltahed  tight,  bounds  for  the  quotient  of  two  difforont 
lull -noma  of  the  sumo  tri -diagonal  matrix  J  ,  the  author  observes  that 
these  bound  a  could  bo  of  uuti  in  un  on'oi’-nna  lyiiJ  s  prov  ided  a  suitable 
algorithm  worn  found.  Such  un  algorithm  is  exhibited,  und  its  errors  srn 
thoroughly  accounted  for,  including  t.h#  elToctii  of  uciiling,  over/ under¬ 
flow  and  roundoff.  A  typical  result  in  I, Nit,  on  a  computer  using  rounded 
flouting  point  binary  arithmetic,  the  biggest  a.lgonvuluo  of  J  cun  be 
computed  euully  to  within  2.3  units  In  its  last  place,  and  the  smaller 
eigenvalue!)  will  auffer  absolute  errors  which  are  no  larger.  Those  re¬ 
sults  are  somewhat  stronger  thou  had  boon  known  before. 
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of  the  contractor,  subcontractor,  grantee,  Department  of  Do* 
fenaa  activity  or  othet  organisation  (corporals  author)  Issuing 
lha  report, 

3a,  REPORT  IEGUWTY  CLAIWKICATIONi  Enter  the  over- 
all  security  classification  of  the  report,  Indicate  whether 
"Reetrioted  Dele"  la  Intiludedi  Marking  la  (0  be  In  accord* 
anea  with  appiopitate  security  regulations, 

2ft,  OHQUPi  Automatic  downgrading  la  anaollied  In  DoD  Di> 
restive  8200, 10  and  Armed  Korea*  Industrial  Manual,  Kn(*e 
lha  group  number,  Alsu,  when  applicable,  ahow  lhal  opilona) 
marking*  have  been  ueed  for  tlroup  3  and  Ctroup  4  a»  author* 
laed, 

I,  REPORT  TITLBi  Knter  the  oomplele  report  title  In  all 
capital  lellera.  Tltlae  in  all  uaeea  ahould  be  unolaaalfiedi 
U  a  meaningful  Mile  oannol  be  aelertled  without  olaeelfloa* 
tlon,  ahow  title  olaaalfluatlon  In  all  uapltala  in  parenthaala 
Immediately  following  the  title, 

4,  DKICklPTIVK  NOTKHi  If  appropriate,  enter  the  type  of 
repurt,  a,  a,,  Interim,  prograaa,  summary,  annual,  or  final, 
tllva  the  inoluaive  dates  when  a  specific  reportini  period  la 
covered. 

It,  AUTIK)H(0)i  Knter  the  name(a)  of  authoKa)  aa  ahown  on 
or  In  the  ration,  Bnltt  teat  name,  flrat  name,  middle  Initial, 

If  military,  ahow  rank  and  branch  of  aerviue,  The  name  of 
the  principal  author  la  an  abaolute  minimum  requirement, 

ft,  HKl'OHT  DATlb  Knter  the  data  of  the  report  ae  day, 
month,  yean  or  month,  yaar,  If  mora  Ilian  one  dale  eppeere 
on  llte  report,  uee  date  of  publication, 

7a,  TOTAL  NUMIIKK  OK  PAOKHi  The  total  page  count 
ahould  follow  normal  pagination  prooedurea,  he,,  enter  the 
number  of  pegee  uonlebting  Information, 

7 ft,  NUMIIKK  OK  NKKKNKNCK*  Unler  the  total  number  of 
referenueu  uiled  In  the  report, 

tta,  CONTRACT  OK  UNANT  NUMIIKK, 1  If  appropriate,  enter 
the  applloabla  number  of  the  contract  or  grant  under  which 
(he  report  wee  written 

16,  go,  tt  «t f,  PKOJlfCT  NUMDKKi  Knter  the  appropriate 
military  department  identlfiuation,  auoh  aa  project  number, 
■ubprojetil  numbai,  ayatem  number*,  leak  number,  etc, 

Da,  ORIOINATOM'N  MICPOKT  NUMUUR(B)i  Knter  the  offl- 
uf «1  rettort  numbar  by  whiuh  the  document  will  be  Identified 
end  controlled  by  the  originating  nativity,  Tlda  number  muet 
be  unique  to  thia  report, 

Oft,  OT I  IKK  REPORT  NlJhUilCW(t*)l  If  the  report  ha*  been 
aaalgned  any  other  report  number*  (allhar  by  lha  orlilnalor 
or  by  lha  a/wnaor),  alao  enter  thl*  number(a), 

10,  AVAIL AU1LITY/LIM1TATION  NOTICKKi  Knter  any  lim* 
ttatlona  on  further  dleeemlnatlon  of  the  report,  other  than  thoa# 


-ItnpSlUf  bV  security  classification.  uatng  standard  atatemtnta 
auch  aai 

(1)  “Qualified  requester*  may  obtain  eopUa  of  this 
report  from  DUG" 

(2)  "Kortdgn  announeemenl  and  dleeemlnatlon  of  thle 
report  by  DUG  le  net  authorised," 

(&)  "U,  g,  Oovernment  agent) lee  may  obtain  ooplea  of 
thla  repoit  directly  from  DUG,  Other  qualified  DDG 
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thla  repoit  directly  I  row  DUG 
users  ahull  requeal  through 


"U,  M,  military  agenulea  may  obtain  ooplea  of  this 
report  directly  from  DDG  Other  queltliml  ueere 
•hall  requeil  through 


(ft)  "All  distribution  of  Dili  report  la  controlled,  Qunl* 
tfled  DDG  ueere  shall  requeal  through 

If  <h*  taport  has  haan  furnished  lo  the  Office  of  Technical 
darvlcea,  Department  of  Gcmmerce,  for  sale  io  lha  public,  Indl* 
cats  this  fact  and  enter  lha  price,  if  known, 

11,  KUKPLKMKNTAHY  NOTKii  Use  for  additional  *  apt  ana* 
lory  nolae, 

12,  IKONkOMlNlt  MILITARY  AOTIVITYl  Urder  the  name  of 
the  departmental  project  office  or  lebnretory  aponaurlng  (pay¬ 
ing  for)  (ha  research  and  development,  Include  address, 

13,  AllHTRAGTi  Knter  an  abelreul  giving  e  brief  end  factual 
summary  of  the  douument  Indicative  of  (he  report,  even  though 
it  may  alao  appear  eieewhere  In  the  body  ol  (he  technical  re* 
noil,  If  additional  apace  le  required,  a  continuation  sheet  ahall 
be  attached. 

It  la  highly  desirable  that  the  abstract  of  olaaalfled  reports 
be  unolaaelned,  Ktult  paragraph  of  (he  abstract  shall  end  with 
an  Indication  of  th*  military  aeuurity  classification  of  the  In* 
formation  in  tha  paragraph,  represented  as  < ti),  fa),  fC),  at  ((/>, 

’liters  la  no  limitation  on  the  length  of  the  abstract,  How* 
ever,  tha  suggested  length  la  from  ISO  to  228  woide, 

14,  KUY  WOHD81  Key  words  era  technically  meaningful  terms 
ur  short  phrase*  Hurt  chareuterlse  a  report  end  may  b*  used  aa 
iitdaa  entries  for  uaialoglcg  the  report,  Key  word*  muet  bo 
selected  eo  the!  no  security  claeaifluetlon  is  required,  tdentl* 
fiere,  such  sa  equipment  mcdel  iteulgnation,  trad*  name,  military 
project  coda  name,  geographic  location,  may  Ire  used  aa  bey 
words  hut  wilt  b#  followed  by  an  indluaUon  of  tauhnica!  con¬ 
test,  The  aaelgnmanl  of  links,  relaa,  and  waighta  is  optional, 
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